1 Macula Densa Project

1.1 Objectives

Three Main Goals of this File
Produce Cleaner looking code.
Identify the amount of clusters there are
Identify the top genes expressed in each of the clusters

2 Loading in Data sets + Library packages.

## 
## The downloaded binary packages are in
##  /var/folders/7n/x74qctp91rng390gx0z9hmd80000gn/T//RtmpIXQzec/downloaded_packages
options(future.globals.maxSize = 74 * 1024^3) # 55 GB
getOption("future.globals.maxSize") #59055800320
## [1] 79456894976
load(here("jk_code", "JK_cleanMD_NEW.rds"))

3 Analyzing DATASET

SO4@meta.data$sample <- factor(
  SO4@meta.data$sample, 
  levels = c("SO1","SO4","SO2","SO3"))
  
  
DimPlot(SO4,split.by ="sample",group.by = "sample")

DimPlot(SO4,group.by = "subclass2_MD",split.by = "sample")

DimPlot(SO4,group.by = "subclass_MD",split.by = "treatment")

DimPlot(SO4,group.by = "subclass2_MD",split.by = "treatment")

DimPlot(SO4)

3.1 Viewing for different DEGs

markers.to.plot1 <- c(
  "Atf3",     # 
  "Egr1",     # 
  "Fos",      # 
  "Jun",      #
  "Junb",     #
  "Pappa2",   #
  "Cxcl10",   # 
  "Cldn19",   # 
  "Krt7",     #
  "Egf",      #
  "Ptger3",
  "Ckb",
  "Mcub",
  "Fabp3",
  "Foxq1",
  "Vash2",
  "Pamr1",
  "Vegfa"
)

                      
DotPlot(SO4,
features = markers.to.plot1,
dot.scale = 8,
dot.min = 0,
scale = FALSE,
scale.max = 100,
scale.min = 0,
col.min = -2.5,
group.by = "seurat_clusters",
col.max = 2.5)+
coord_flip()

DimPlot(SO4,group.by = "seurat_clusters")

DimPlot(SO4,split.by = "treatment",group.by = "seurat_clusters")

library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:AnnotationDbi':
## 
##     select
## The following object is masked from 'package:IRanges':
## 
##     slice
## The following object is masked from 'package:S4Vectors':
## 
##     rename
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
dim_plot <- DimPlot(SO4,split.by = "treatment",group.by = "seurat_clusters") +
  ggtitle("UMAP Dimensional Reduction")
interactive_plot <- ggplotly(dim_plot)

# Display the interactive plot
interactive_plot

3.2 Finding Markers for each MD Subtype

Idents(SO4) <- "subclass2_MD"

MD_DEGs <- FindAllMarkers(SO4, 
                          only.pos = TRUE,
                          logfc.threshold = 0.1,
                          min.pct = 0.1,
                          min.diff.pct = 0.1,
                          return.thresh = 0.05)
## Calculating cluster type_1
## Warning: The `slot` argument of `GetAssayData()` is deprecated as of SeuratObject 5.0.0.
## ℹ Please use the `layer` argument instead.
## ℹ The deprecated feature was likely used in the Seurat package.
##   Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: `PackageCheck()` was deprecated in SeuratObject 5.0.0.
## ℹ Please use `rlang::check_installed()` instead.
## ℹ The deprecated feature was likely used in the Seurat package.
##   Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## For a (much!) faster implementation of the Wilcoxon Rank Sum Test,
## (default method for FindMarkers) please install the presto package
## --------------------------------------------
## install.packages('devtools')
## devtools::install_github('immunogenomics/presto')
## --------------------------------------------
## After installation of presto, Seurat will automatically use the more 
## efficient implementation (no further action necessary).
## This message will be shown once per session
## Calculating cluster type_2
## Calculating cluster type_3
## Calculating cluster type_4
MD_DEGs <- MD_DEGs %>% arrange(desc(avg_log2FC))

write.xlsx(MD_DEGs, file = here("jk_code", "DEGs_MD_SUBclass.xlsx"), rowNames = TRUE)

# save as excel file

MD_DEGs %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1, p_val_adj < 0.05) %>%
    slice_head(n = 10) %>%
    ungroup() -> top10


# save as excel file

DoHeatmap(SO4, features = top10$gene) + NoLegend()
## Warning in DoHeatmap(SO4, features = top10$gene): The following features were
## omitted as they were not found in the scale.data slot for the SCT assay: Oasl2

top10

3.3 Markers for each MD Sub-SUb type

Idents(SO4) <- "subclass_MD"

MD_DEGs2 <- FindAllMarkers(SO4, 
                          only.pos = TRUE,
                          logfc.threshold = 0.1,
                          min.pct = 0.1,
                          min.diff.pct = 0.1,
                          return.thresh = 0.05)
## Calculating cluster type_1
## Calculating cluster type_2a
## Calculating cluster type_2b
## Calculating cluster type_3
## Calculating cluster type_4
MD_DEGs2 <- MD_DEGs2 %>% arrange(desc(avg_log2FC))

write.xlsx(MD_DEGs2, file = here("jk_code", "DEGs_MD_SUBSUBclass.xlsx"), rowNames = TRUE)

# save as excel file

MD_DEGs2 %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1, p_val_adj < 0.05) %>%
    slice_head(n = 10) %>%
    ungroup() -> top10


# save as excel file

DoHeatmap(SO4, features = top10$gene) + NoLegend()
## Warning in DoHeatmap(SO4, features = top10$gene): The following features were
## omitted as they were not found in the scale.data slot for the SCT assay: Oasl2,
## Atp6v0d2, Irx3

3.4 Finding markers between each type now

3.4.1 Type 1a Markers

type_1 <- c("Pappa2","Ramp3","Itga4","Aard")

VlnPlot(SO4,type_1,)
## Warning: The `slot` argument of `FetchData()` is deprecated as of SeuratObject 5.0.0.
## ℹ Please use the `layer` argument instead.
## ℹ The deprecated feature was likely used in the Seurat package.
##   Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

3.4.2 Type 2a Markers

type_2a <- c("Egf","Umod","Cxcl12","Ckb")

VlnPlot(SO4,type_2a)

3.4.3 Type 2b Markers

type_2b <- c("Foxq1","Slc9a3","Sult1d1","Tmem52b")

VlnPlot(SO4,type_2b)

3.4.4 Type 3 Markers

type_3 <- c("Jun","Fos","Egr1","Hspa1a")

VlnPlot(SO4,type_3)

3.4.5 Type 4 Markers

type_4 <- c("Cxcl10","Isg15","Ifit1")

VlnPlot(SO4,type_4)

3.5 Top Genes for each type

type_markers <- c(type_1, type_2a, type_2b, type_3, type_4)

                      
                      
DotPlot(SO4,
features = type_markers,
dot.scale = 8,
dot.min = 0,
scale = FALSE,
scale.max = 100,
scale.min = 0,
col.min = -2.5,
col.max = 2.5)+
coord_flip()

4 Save Object

save(SO4, file = here("jk_code", "SO4_SUB_NEW_analysis.rds"))